!***********************************************************************
!*   Программа, вычисляющая методом сеток характеристики надкритических 
!*   движений для соответствующих чисел Прандтля (PR) и Рэлея (RA) 
!*   в ячейке Хеле - Шоу.
!***********************************************************************
!*     DT      - Шаг по времени
!*     TIME    - Текущее время
!*     L       - Длина ячейки
!*     H       - Высота ячейки  
!*     HX      - Шаг по координате X
!*     HY      - Шаг по координате Y
!*     Q(I,J)  - Поле вихря
!*     F(I,J)  - Функция тока       
!*     T(I,J)  - Поле температуры
!*     Q1(I,J) - Дополнительный массив для вихря
!*     F1(I,J) - Дополнительный массив для функции тока
!*     T1(I,J) - Дополнительный массив для температуры 
!*     FMAX    - Максимальное значение функции тока в массиве
!*     FMIN    - Минимальное значение функции тока в массиве
!*     NU      - Число Нуссельта
!***********************************************************************
      PROGRAM NET
      IMPLICIT REAL*8 ( A - Z )
      INTEGER I, J, M, FLAG, N1, N2, M1, M2, K
      PARAMETER ( M1 = 151, M2 = 79 )
      DIMENSION Q(M1,M2),  T(M1,M2), F(M1,M2)
      DIMENSION Q1(M1,M2),  T1(M1,M2), F1(M1,M2)
      LOGICAL EX
      COMMON /BLK/ DT, HX, HY, N1, N2, EP

      RA=11.0d0      !   RA   -  число Рэлея
      L=45.0d0       !   L    -  длина ячейки
      H=20.0d0       !   H    -  высота ячейки
      PR=20.0d0      !   PR   -  число Прандтля
      EP=0.00001d0   !   EP   -  погрешность для метода ПВР
      EN=0.0d0       !   EN   -  погрешность выхода на стационар
      FLAG=2         !   FLAG -  временная зав-ть функции тока
      TEXIT=20000    !   TEXIT - время выхода
      N1=151         !   N1   -  число узлов вдоль x
      N2=79          !   N2   -  число узлов вдоль у

5     FORMAT ( A1 )
      CLOSE ( UNIT = 5, STATUS = 'KEEP' )   
!******   Максимум функции тока в зависимости от времени
      IF (FLAG.EQ.0) THEN
        open (unit = 1000, file = 'psi1.dat', status = 'new')
        open (unit = 2000, file = 'psi2.dat', status = 'new')
        open (unit = 3000, file = 'psi3.dat', status = 'new')
        open (unit = 4000, file = 'psi4.dat', status = 'new')
        open (unit = 5000, file = 'psi5.dat', status = 'new')
      ENDIF
!******   Вычисление вспомогательных величин
      PI = 3.14159265358979
      HX = L/DBLE( N1-1 )
      HY = H/DBLE( N2-1 )
      HX2 = HX*HX
      HY2 = HY*HY
      DT = MIN(HX2,HY2)/5.6
!******   Присвоение начальных значений массивам
      DO I = 2, N1 - 1
        DO J = 2, N2 - 1 
          Q(I,J) = 0.0d0
          T(I,J) = 0.0d0
          F(I,J) = 0.0d0
        ENDDO
      ENDDO
!******   Задание граничных условий
      DO I = 1, N1
         T(I,1) = 0.0d0
         T(I,N2) = 0.0d0
         F(I,1) = 0.0d0
         F(I,N2) = 0.0d0
         Q(I,1) = 0.0d0
         Q(I,N2) = 0.0d0
      ENDDO
      DO J = 2, N2 - 1
         T(1,J) = T(2,J)
         T(N1,J) = T(N1-1,J)
         F(1,J) = 0.0d0
         F(N1,J) = 0.0d0
         Q(1,J) = 0.0d0
         Q(N1,J) = 0.0d0  
      ENDDO
!******   Дополнительные граничные условия
      DO I = 1, N1
         T1(I,1) = 0.0d0
         T1(I,N2) = 0.0d0
      ENDDO		    
!******   Возмущение полей
         T(9,5) = -1.9d0
         T(71,5) = -1.9d0
         F(9,5) = -1.5d0
         F(71,5) = 1.5d0
!******   Вычисление коэффициентов
        A1 = 2.0*DT/(3.0*PI*PR*HX*HY)
        A2 = DT/HX2
        A3 = DT/HY2
        A4 = 1 - 0.25*DT*PI*PI
        A5 = 2.0*RA*DT/(PI*HX)
        A9 = 0.5*DT/(PI*PR*HX*HY)
        A10 = DT/(PI*PR*HX)
        A11 = A2/PR
        A12 = A3/PR
!******   Начальное значение для максимума функции тока
      K = 0
      FP = 0.0d0
      PRINT 25, RA, PR, EP
25    FORMAT (' RA=', 1P1E11.4,' PR=', 1P1E9.2,' EP=', 1P1E13.6 )
      
      DO M = 1, 90000000
         TIME = DT*M
         
C$OMP PARALLEL
C$OMP DO PRIVATE(J, S1, S2, D1, D2, D3, D4, D5, D6, D7, D8)
         DO I = 2, N1 - 1
            DO J = 2, N2 - 1

              S1 = F(I+1,J) - F(I-1,J) 
              S2 = F(I,J+1) - F(I,J-1)
              D1 = Q(I,J+1) - Q(I,J-1)
              D2 = Q(I+1,J) - Q(I-1,J)
              D3 = Q(I+1,J) - 2.0*Q(I,J) + Q(I-1,J)
              D4 = Q(I,J+1) - 2.0*Q(I,J) + Q(I,J-1)
              D5 = T(I,J+1) - T(I,J-1)
              D6 = T(I+1,J) - T(I-1,J)
              D7 = T(I+1,J) - 2.0*T(I,J) + T(I-1,J)
              D8 = T(I,J+1) - 2.0*T(I,J) + T(I,J-1)

              Q1(I,J) = A1*( S1*D1 - S2*D2 ) + A2*D3 + A3*D4 + A4*Q(I,J) - A5*D6 
              T1(I,J) = T(I,J) + A9*( S1*D5 - S2*D6 ) - A10*S1 + A11*D7 + A12*D8
            ENDDO
         ENDDO
C$OMP END DO 
C$OMP END PARALLEL
         
!******   Вычисление новых условий на боковых границах
        DO J = 2, N2 - 1
           T1(1,J) = T1(2,J)
           T1(N1,J) = T1(N1-1,J)
        ENDDO
!******   Функции тока скорости и пульсационной компон. скорости
        CALL PVR ( F, F1, Q1 )
        write (1000, 300) time, F(18, 34)
        write (2000, 300) time, F(133, 34)
        write (3000, 300) time, F(44, 56)
        write (4000, 300) time, F(107, 56)
        write (5000, 300) time, F(76, 56)
        DO I = 2, N1 - 1
           DO J = 2, N2 - 1 
              F(I,J) = F1(I,J)
           ENDDO
        ENDDO
        FMAX = 0.0d0
        DO I = 2, N1 - 1
          DO J = 2, N2 - 1
             FMAX = MAX( F(I,J), FMAX )
          ENDDO
        ENDDO
        FMIN = FMAX
        DO I = 2, N1 - 1
          DO J = 2, N2 - 1
             FMIN = MIN( F(I,J), FMIN )
          ENDDO
        ENDDO
        EXTR = MAX( ABS( FMIN ), ABS( FMAX ) )
!******   Переобозначение массивов
        DO I = 2, N1 - 1
           DO J = 2, N2 - 1
              Q(I,J) = Q1(I,J)
              T(I,J) = T1(I,J)
           ENDDO
        ENDDO
        DO J = 2, N2 - 1
           T(1,J) = T1(1,J)
           T(N1,J) = T1(N1,J)
        ENDDO  
!******   Вычисление вихря на границах слоя по формуле Тома
        DO I = 1, N1
          Q(I,1) = 2.0*F(I,2)/HY2
          Q(I,N2) = 2.0*F(I,N2-1)/HY2
        ENDDO
        DO J = 1, N2
          Q(1,J) = 2.0*F(2,J)/HX2
          Q(N1,J) = 2.0*F(N1-1,J)/HX2
        ENDDO
!******   Нахождение погрешности
        FN = ABS( EXTR - FP )
        IF ( EN.GT.FN ) EXIT
        IF ( TIME.GT.TEXIT ) EXIT ! Ручной выход
        FP = EXTR
!******   Запись временной зависимости максимума функции тока
        IF ( FLAG.EQ.0 ) THEN
          K = K + 1
          IF ( K.EQ.6 ) THEN
!******   Вычисление числа Нуссельта
          NU = 0.0d0
          DO I = 1, N1
             NU = ( T(I,1) - T(I,2) )*HX/HY + NU
          ENDDO
             NU = 2.0*NU/(PI*L)
             PRINT *, TIME, EXTR
             WRITE ( 7, 300 ) TIME, NU
               D = T(8,20) + T(8,30) - T(22,20) - T(22,30)
             WRITE ( 13, 300 ) TIME, EXTR, D
300          FORMAT( 1PE13.6,'  ', 1PE13.6,'  ', 1PE13.6 )
          K = 0
          ENDIF
        ENDIF
        
      ENDDO
    
      CLOSE ( UNIT = 7, STATUS = 'KEEP' )
      CLOSE ( UNIT = 13, STATUS = 'KEEP' )
      IF ( FLAG.NE.1 ) THEN
!******   Запись в файл поля функции тока
      OPEN ( UNIT = 17, FILE = 'FUN.DAT', STATUS = 'NEW' )
      DO I = 1, N1
        DO J = 1, N2
          WRITE ( 17, 640 ) I, J, F(I,J)
640       FORMAT( I3,', ', I3,', ', 1PE13.6 )
        ENDDO
      ENDDO
      CLOSE ( UNIT = 17, STATUS = 'KEEP' )
!******   Запись в файл поля температуры
      OPEN ( UNIT = 27, FILE = 'TEM.DAT', STATUS = 'NEW' )
      DO I = 1, N1
        DO J = 1, N2
          WRITE ( 27, 750 ) I, J, T(I,J)
750       FORMAT( I3,', ', I3,', ', 1PE13.6 )
        ENDDO
      ENDDO
      CLOSE ( UNIT = 27, STATUS = 'KEEP' )
      ENDIF

      STOP
      END